Scenario

Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.

Load packages

library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement 
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed 
library(bayesplot) # needed for MCMC diagnostics 
library(DHARMa) # model validation 
library(ggdist) # partial plots 
library(tidybayes) # partial plots 
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans 
library(ggpubr) # arranging figures 
library(mgcv) # GAM

Set working directory

knitr::opts_knit$set(root.dir=working.dir)

Import data

This data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty

lat_resp_dat <- read_delim("C:/Users/jc527762/OneDrive - James Cook University/PhD dissertation/Data/Chapter4_Latitude/import_files/lat_resp_dat_no_outliers.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(area = col_skip(), rate.a.spec = col_skip()), 
    trim_ws = TRUE)

Data manipulation

rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |> 
  mutate(dev.temp = as.factor(dev.temp), 
         replicate = str_sub(sampleID, -1,-1), 
         population = factor(population)) |>
                                           # number of observations = 5758
  filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
         sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
         sampleID != "60.LCKM.152.30.1"  # 5637 - 76 = 5542
  ) |> 
  filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes) 
  group_by(sampleID) |> 
  mutate(max_value_index = which.max(rate.output), 
         row_number = row_number(), 
         max.rate = max(rate.output), 
         low.rate = mean(rate.output[order(rate.output)[1:10]]), 
         net.rate = max.rate - low.rate) |> 
  filter(row_number <= max_value_index) |>
  ungroup() |> 
  mutate(region = factor(region), 
         dev.temp =factor(dev.temp), 
         level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
                           tank >= 200 & tank <= 299 ~ 2,
                           tank >= 300 & tank <= 399 ~ 3,
                           TRUE ~ NA_real_)), 
         female = factor(female))

Exploratory data analysis

2nd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

3rd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

Fit model [random factors]

First we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.

Hypothesis test will include:

  1. Null model
  2. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK)
  3. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL)
  4. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)
  5. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)+ (1| CLUTCH_ORDER)
lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output))) 
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") +  
  prior(student_t(3, 0, 195), class = "sigma")

Models

Null

f.model.null <- bf(rate.output ~ 1, 
                   family = gaussian()) 

model.null <- brm(f.model.null,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model.null, file="./03.CTmax_resp_data_analsyis_files/model.null.RDS")

Model1

f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank), 
                          family=gaussian())

model1 <- brm(f.model1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1, file="./03.CTmax_resp_data_analsyis_files/model1.RDS")

Model2

f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level), 
                          family=gaussian())

model2 <- brm(f.model2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model2, file="./03.CTmax_resp_data_analsyis_files/model2.RDS")

Model3

f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population), 
               family=gaussian())

model3 <- brm(f.model3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model3, file="./03.CTmax_resp_data_analsyis_files/model3.RDS")

Model4

f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order), 
                        family=gaussian())

model4 <- brm(f.model4,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model4, file="./03.CTmax_resp_data_analsyis_files/model4.RDS")

LOO

LOO(model.null,model1, model2,model3,model4) 
## Output of model 'model.null':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -31469.1  78.1
## p_loo         2.5   0.4
## looic     62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.6  99.6
## p_loo        57.6   3.1
## looic     61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.5
## p_loo        57.6   3.0
## looic     61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.7  99.6
## p_loo        58.1   3.2
## looic     61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model4':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30662.5  99.4
## p_loo        57.5   3.0
## looic     61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##            elpd_diff se_diff
## model4        0.0       0.0 
## model1       -0.1       0.4 
## model2       -0.2       0.4 
## model3       -0.2       0.4 
## model.null -806.6      48.2

Fit model [fixed factors - polynomials]

Prior investigation

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output))) 
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors <- prior(normal(430, 0.25), class="Intercept") + 
  prior(normal(0, 40), class="b")
  prior(student_t(3, 0, 195), class = "sigma")

Models

Model1.p1

Model1.p1

f.model1.p1 <- bf(rate.output ~ 1 + 
                         dev.temp*region + t +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p1 <- brm(f.model1.p1,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p1, file="./03.CTmax_resp_data_analsyis_files/model1.p1.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.987, p-value = 0.656
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p1$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p1$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p1$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p1$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p1$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

model1.p2

model1.p2

f.model1.p2 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(t, 2) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p2 <- brm(f.model1.p2,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p2, file="./03.CTmax_resp_data_analsyis_files/model1.p2.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98869, p-value = 0.744
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p2$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p2$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p2$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p2$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p2$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

model1.p3

model1.p3

f.model1.p3 <- bf(rate.output ~ 1 + 
                         dev.temp*region + poly(t, 3) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p3 <- brm(f.model1.p3,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p3, file="./03.CTmax_resp_data_analsyis_files/model1.p3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98619, p-value = 0.576
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p3$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p3$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p3$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p3$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p3$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

LOO

LOO(model1.p1, model1.p2,model1.p3) 
## Output of model 'model1.p1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -29955.1  94.0
## p_loo        60.1   3.1
## looic     59910.2 188.0
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30526.2 100.6
## p_loo        58.4   3.1
## looic     61052.4 201.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -30519.9 100.2
## p_loo        59.1   3.2
## looic     61039.7 200.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##           elpd_diff se_diff
## model1.p1    0.0       0.0 
## model1.p3 -564.8      31.2 
## model1.p2 -571.1      31.4

Fit model [fixed factors - GAM]

GAM model

f.model1.gam <- bf(rate.output ~ 1 + 
                         dev.temp*region + s(t) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.gam <- brm(f.model1.gam,
              data = lat_resp_dat2, 
              prior = rate.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.gam, file="./03.CTmax_resp_data_analsyis_files/model1.gam.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gam, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gam, ndraws=250, summary=FALSE)
model1.gam_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gam_resids) ; testDispersion(model1.gam_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.98442, p-value = 0.6
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gam$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gam$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gam$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gam$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gam$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.